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Abstract: Mode properties of circularly symmetric waveguides with one 
special type of bianisotropy are studied using finite element approach. We 
find that the polarization degeneracy in circularly symmetric waveguides 
can be eliminated, by introducing intrinsic crossing coupling between 
electric and magnetic moments in the constituent units of the waveguide 
media. Breaking the polarization degeneracy in high order mode groups is 
also confirmed numerically. With the bianisotropic parameters chosen in 
this work, the x and y-polarized modes remain decoupled. Typically, the 
y-polarized modes remain completely unchanged, while the x-polarized 
modes are turned into leaky modes that are lossy along propagation 
direction. A perturbation model from coupled mode theory is developed 
to explain the results and shows excellent agreement. Such asymmetric 
behavior between different polarizations might be feasible and useful for 
developing compact polarizers in terahertz or mid-infrared regime. 
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1. Introduction 

In metamaterials, the most general description of material properties is bianisotropic consti¬ 
tutive relation, which governs material responses to the electromagnetic (EM) fields, as well 
as the cross terms, i.e., magnetic (electric) dipoles induced by the incident electric (magnetic) 
fields. Conventionally, the cross terms are hidden in the higher order spatial derivatives of the 
electric polarization since there is an ambiguity of splitting the induced current into polarization 
and magnetization in the macroscopic electromagnetics HI . However, by properly selecting the 
‘gauge field’ for the polarization and magnetization, one is able to treat electric and magnetic 
fields on the same footing m, thus be able to employ the duality between electric and mag- 



netic fields to simplify the analysis significantly. It was found that most of moving media are 
bianisotropic by Kong W, who invented the concept of bianisotropy to describe the electro¬ 
magnetic properties of such kind of materials. 

One might expect exotic behavior of EM waves, or novel devices may exist iiaibi, due to 
the extra degree of freedom from materials, i.e., the EM coupling in constituent units or the 
building blocks. Indeed, extensive research activities, including experimental efforts CHi and 
theories UMIl, have been carried out in conceiving and synthesizing bianisotropic metama¬ 
terials. As a typical sub-class of bianisotropic media, chiral media are reciprocal, and can be 
useful in circular dichroism spectroscopy QSl, and polarization control. On one hand, there are 
considerable efforts on the fundamental level to generate chiral light, or super-chiral light from 
chiral or achiral photonic structures, i.e., helix metallic structures HKTh) and planar photonic 
lattices lEiini. On the other hand, Engheta and others seek applications using chiral effects, 

1. e., chirowaveguides II18II19II and chiral fibers ll20l . Chirowaveguides proposed by Engheta are 
realized by filling cylindrical waveguides with isotropic chiral media, while chiral fibers are 
experimentally feasible with current fiber technology, i.e., by microforming glasses or deco¬ 
rating glass surfaces to create helical structures. Inspired by Engheta’s pioneering work 03, 
chirowaveguides are examined by many authors imHu. We note that in all the aforementioned 
work, the chiral media are isotropic. As for light, it corresponds to a macroscopic model where 
chiral molecules are randomly distributed in the host, while for microwave the counterpart is 
the wire helices with random orientations in the host media. However, most of the available 
experiments in the field of chiral metamaterials HKTh) suggest that it is very likely that the 
synthesized chiral media have strong bianisotropy. To the best of our knowledge, an efficient 
tool of analyzing bianisotropic waveguides is in lack in literature. 

In this paper, we extend the aforementioned work on chirowaveguides to the regime where 
the media can be bianisotropic using a finite element approach. Particularly, we study the mode 
properties of bianisotropic waveguides with a special emphasis on how the guided x- and y- 
polarized modes are impacted by bianisotropy. In contrast to Engheta’s findings that TE and 
TM modes are entangled together, we found that our bianisotropic waveguides essentially func¬ 
tion as linear polarizers. The elimination of polarization degeneracy ll25l typically refers to two 
effects, i.e., birefringence and different propagation lengths of the x- and y-polarized modes. 
The difference between real parts of effective modal indices (Re{neff)) accounts for the bire¬ 
fringence, while different lm{nf,ff) leads to different propagation lengths. 

The paper is organized as follows. In Section 2, we outline the geometry under investiga¬ 
tion and give the theoretical foundation of how the finite element approach can be applied for 
bianisotropic waveguides. In Section 3, we study the mode dispersion and give a coupled mode 
model to explain our results. Einally, the paper is concluded in Section 4. 

2. Bianisotropic waveguides and theoretical foundation of FEM modelling 

2.7. Bianisotropic waveguides 

The bianisotropic waveguide studied in this work is sketched in Eig. [Ha). The constitutive 
relation can be given as follows, 

D=-eE + x,i,H, 

B = ilH + XheE, 

where £ = e^Er, p. = liop,, Xhe = VeoMoZL’ Xeh = Reciprocity imposes the con¬ 

straints on the material parameters as given by £ = £^, p = p^, Xeh — ~Xj,e- lossless 
bianisotropic media, Xeh 1® purely imaginary. Bianisotropic media in principle can be realized 
by aligning electrically small magnetoelectric inclusions along particular directions. EigurelHb) 
andfljc) show two examples of aligned dipoles from small structures, i.e. split-ring resonators 
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Fig. 1. Schematic diagram of a circularly symmetric waveguide (a) which contains bian- 
isotropic medium in the core layer by including structures with split ring resonator (SRR) 
arrays introducing ;(i 2 (b), and structures with helix arrays introducing (c). 


(SRRs) and helix arrays, respectively, filling a typical block (dashed region) shown in Fig. [Ha). 
It should be noted that effective constitutive parameters of SRRs and helix arrays depend on 
many factors such as geometric shape, lattice constant, and incident polarization Il26l - l28l . The 
bianisotropy could be very dispersive as well. 

In our case, we study bianisotropic waveguides operating at a single frequency. The dimen¬ 
sions of SRRs or helix arrays should be far smaller than the wavelength, which implies that 
structure details, i.e., geometric shapes and lattice constant of the meta-atoms, are averaged out, 
except the bianisotropic response built from the aligned meta-atoms. Under this assumption, 
applying a magnetic field along the y-direction on split ring resonators (SRRs) will induce elec¬ 
trical currents flowing over the SRRs along clock or anti-clock direction, which subsequently 
generate electric dipoles orientated along .r-direction, as shown in Fig.ITJb). This corresponds 
to an off-diagonal term Xi 2 in Xe/t- Figure Ec) shows a typical structure that can be used to 
generate the diagonal term jn in 

2.2. Theoretical foundation ofFEM modelling 

We describe the theoretical foundation of the finite element method that is employed in this 
work. We take a time harmonic dependence for the EM waves throughout this paper. One 
can rewrite Eq. (E and insert it into source-free Maxwell’s equations, 

V X £ = -i(o[fLX-jD + {X,!,, - ilXfh£)E], .21 

^x{x;i!D-x;i!£E) = icoD. 

Combining the two equations in Eq. E, one can calculate the eigenmodes of the bianisotropic 
waveguides as given in the following equation, 

-V X jU^'V X E{r,co) +koCi'V xE{r,co)+k(f7 x C 2 E{r,co) +kQ{di +d 2 )E{r,co) = 0, (3) 

where di = £r, d 2 = -Xehf^^'Xhe’ b 2 = -iflf^Xhe’ and ko = coo^eoHo- Due to 

the translation symmetry of the waveguides, we can expand the modes in combination of plane 

waves along z-direction, i.e., E{r) = J E(x,y,l5)e^’^^dp, where £(x,y,j3) = xEx{x,y,p) + 

yEy{x,y,j5) + zE^{x,y,j5). Substituting E{r) into Eq. E in conjunction with the test function 
of F{x,y^z) = F(x,y,j 32 )e^'^ 2 z^ vve can get the variational form of Maxwell’s equation for the 


















waveguide modes give by 2-dimensional (2D) integration over the transverse plane (S), 

L = — JJdxdy—CurlE{x,y,j5) ■CurlF*{x,y,P) + F*{x,y,l5) ■ [—n x CurlE{x,y,p)]dl 

S dS 

+ jj dxdyh)[c\^ ■KC2 + h){d\+d 2 )\E{x,y,^)-F*{x,y,P), 
s 

(4) 

where CurlE{x,y,li) = [x{^ +y{-^ - iliE^)+z{^ - ^)]. For guided 

modes, the fields in the transverse plane decay to zero on the boundaries, as the modelling 
domain is relatively large, i.e., several tens of wavelength. Hence the term of the boundary in¬ 
tegration in Eq. (IHi drops out. As regards to Eq. (|4]i, there are a few facts and implications that 
we intend to discuss. Eirstly, Eq. © is corresponding to the functional of the reduced wave 
equation, and can be implemented in a truncated 2D computation domain, using standard fi¬ 
nite element procedure. Secondly, Eq. @ can be turned into an eigenvalue problem, in which 
the propagation constant j5 and the field solutions £(x,y,j3) can be solved directly. Thirdly, 
Eq. © is capable of handling waveguides consisting of bianisotropic materials, in either core 
or cladding layer, or both of them. In the following section, we will use this numerical tool to 
study a circularly symmetric waveguide, in which the material of the core layer is bianisotropic. 
Eurther details of implementation of EEM calculations can be found in Il29l43^ . 

3. Mode dispersion and coupled mode model 

3.1. Dispersion 

Starting from the simplest case, single mode cylindrical waveguides are considered in the pa¬ 
per. As shown in Eig.[Tla), the geometry is identical to conventional cylindrical waveguide with 
high material index in the core surrounded by air, except that the core layer is bianisotropic. 

/ 0 ;ti2 0 \ 

Explicitly, the cross term is given by x'eh = ~[XheV = i 0 0 0 > which corresponds 

\ 0 0 0 / 

to the scenario sketched in Eig.[Tlb). Rigourously speaking, bianisotropic parameters, permit¬ 
tivity and permeability are not completely independent from each other. Eor instance, Xeh 
homogenous medium as well as that of magnetoelectric point scatters suffer a upper bound Col 
imposed by the magnitude of the electric and magnetic responses, i.e., Xn < In this pa¬ 

per, we approximate the bianisotropic parameters as independent parameters from permittivity 
and permeability, which may hold for the media that are off structure-resonances. Without fur¬ 
ther complications, the permittivity and permeability are approximately taken as constant scalar 
values. 

Eigure|2a) and|2lb) show Ks{iigff) and \m.{ngff) as a function of Xn for three core radii. 
The sizes of the waveguides are chosen so that only fundamental modes with two orthogo¬ 
nal polarizations are supported. It can be seen that y-polarized modes keep almost unchanged, 
while x-polarized modes are altered significantly. In other words, the degeneracy between x- 
and y-polarized modes, a commonly feature of circularly symmetric waveguides, is broken. In 
addition, it is clear that two orthogonally polarized modes are decoupled. Such property is dis¬ 
tinct from isotropic chiral media or chirowaveguides where x- and y-polarization are entangled 
together. Eurthermore, forx-polarized modes, Re(«(.//) drops while \\m{neff) \ increases as Xi 2 
increases. This means that x-polarized modes become lossier as Xi 2 increases. To illustrate this, 
we plot the propagation lengths as a function of Xi 2 when the powers of x-polarized modes are 
attenuated by 1000 times (30dB) compared to the case oflm(ngff) = 0 , as shown in Eig.|2tc). 



1.45 


1.3 

1.15 

1 


(a) 

0.1 

(b) 




X.. 

X 

\ 

\ 

. « . 

'^-0.1 

<u 

c 

^ 03 

\ 



\ 


-0.5 





0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 


X,2 5^12 


100 r 


(C) 


Q. 

O 


80 ^ 

60 

40 

0 

20 » 
0 


0 0.2 0.4 0.6 0.8 


''12 


Fig. 2 . (a) Re.(?igff), (b) lm{ngff), (c) propagation length of x-polarized modes. Symbols 
with (without) lines refers to y- (x-) polarization modes. The color indicates different size 
of the radius, i.e., 0.12Ao (black circles), O.I6A0 (blue triangles) and 0.2Ao (red diamonds) 
respectively, where Ao is vacuum wavelength. We use En = £22 = £33 = 4 , flu = fl22 = 
fl23 = 1 . The other elements in Xg/j is set to be zero. 


Thus high polarization extinction ratio might be realized over a propagation distance of a few 
wavelengths. 

Figure|2]also shows the different behavior of ngff for three different waveguides. In addition, 
x-polarized modes vanish upon certain Xn values. In the following, we develop a coupled mode 
model to explain these phenomena. 

3.2. Coupled mode theory in bianisotropic waveguides 

To understand how the modes are impacted by the bianisotropy, it is useful to examine con¬ 
stitutive relation of the bianisotropic medium firstly. The three components of D and B can be 
written as 

Dx — £o(£ii£x “F Bx = fiofinlix', 

Dy = Eoe22ey, By = fioi—iXnex + iiiiK)', (5) 

Dz = eoEj^ez, Bz = fiDii^'ihz- 

It is clear from Eq. Q that Dy and Bx have no contribution from magnetoelectric coupling in 
our configuration. Therefore, y-polarized modes are not affected since they are dominated by 
Ey and Hx. On the contrary, x-polarized modes are dominated by Ex and Hy, hence they are 
strongly modified due to Xn- In the following, we construct a coupled mode model to study 
how the x-polarized modes evolve as Xi 2 increases. 

To simplify our analysis, we use normalized fields which are given by [e,h.] = [E,ZoH] 

with Zo = Expressing e and h as e{r) = e{x,y,j3)e^‘l^^, h{r) = h{x,y,j5)e^‘l^~, we refor¬ 
mulate the Maxwell’s equation for the bianisotropic waveguide modes as follows 

Vt X - ipz X + xle^^‘‘), (6a) 

V, X - ij5z X + (6b) 

where = e(x,y,j3), = h{x,y,l5), V, = x-^ and j3 is the propagation constant. 

Applying Eq. (|6]l to two modes with different propagation constants and fields, namely, )3o 
and j3 we derive the formula of Aj3 as given by 

-if{exhQ)-ndl + JJ{ikofi,.h)-hl-{koereQ)-e}dxdy //core 

^ “ fjz-iex h*o)dxdy ' ° fjz-{e x hl)dxdy 

(7) 














Fig. 3. (a) Power leakage ratio defined by 7 ] versus defined by Eq. (| 3 l. (b) Real part 

of fij normalized by ex (Cy) for x-{y-) polarized mode versus X\2- Radius of the waveguide 
in (a) and (b) is 0.2Ao. 


where Aj3 = j3 — j3o, t^Xxi = X12 — Xi 2 ,o^ ” the normal vector of the outer boundary. We refer to 
Refs. Il3?l435l for further details on the derivation. It is important to point out the fact that x- and 
y-polarized modes are decoupled, which is in consistency with the dispersion relation shown 
in Fig. I 2 ] and the analysis of Eq. (|5]l. Mathematically, it translates to the fact that there is no 
projection between the fields in Eq. (|7]) from the two modes picked up from two different dis¬ 
persion lines of Eig.|2] This is in marked contrast to the chirowaveguides discussed by Engheta. 
In the decoupled case, Eq. © is exact and shows quantitatively how p evolves as Xn varies. 
As a trivial application of Eq. © to y-polarized modes, we immediately find that their prop¬ 
agation constants keep unchanged, since Xhe (Xeh^ picks up Ex (Hy) components, which 
are far smaller than Ey (Hx). The derivation of the coupled mode theory in coupled cases are 
similar. However, the final form of Eq. © shall be modified in a matrix form, whose elements 
are determined by the projection between different polarization modes. 

3.3. Power leakage of x-polarized modes 

We proceed to discuss how the x-polarized modes are impacted by the variation of Xn- In case 
of no variation, i.e., Axn = 0, and h ^ ho, e = cq, Aj3 = 0, Eq. © is reduced to the following 
relation given by 

X ■ndl = jj {(koA,0 • Kr - ikoSriel^r) ■ ef}dxdy. (8) 

It is well known that for properly propagating waveguide modes, namely, perfectly guided 
modes, the energy stored in electric and magnetic fields shall be balanced. As evident from 
Eq. ©, the difference of the stored energy between electric and magnetic part equals to the 
leakage of power flux in the transverse plane. Indeed, Eig. |3a) shows the leakage ratio T] as a 
function of Xn where rj is defined as 

If {(koAr^o) ■ - (kpe^gg) ■ epjdxdy 
ff {(kofixf>o) ■ fiQ + (ko£reg) ■ eojdxdy' 

It can be seen that t] increases (t] « 0 ) with Xn for x (y)-polarized modes, which is consistent 
with changes of Im(«e//) of x (y)-polarized mode shown in Eig. |2l 

The power leakage of the x-polarized modes can be understood by analyzing the field compo¬ 
nents, typically z components, in a simple yet instructive manner. We can find e^ from transverse 
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components of electric and magnetic fields, i.e., 
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= 2x2, 
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the real part of which increases as Im(j3) gets large. Figure [^b) shows the ratio of for 

x-polarized mode {ex dominated in this case) and | \ for y-polarized mode (e, dominated) 

at the center of the core layer as a function of Xn- For y-polarized mode the ratio remains 
constantly small, while it increases continuously with Xn for x-polarized mode. We note the 
hj component bears the same trend as e^, which is not shown here. The increase in Re(ej.) 
gives rise to the power leakage in the transverse plane, as evident from the non-zero transverse 
components of the real part of Poynting vector given by 


Px = Re{/[evlm(/zj) + /!vlm(e,)] -f [eyRe{h,) - /r,.Re(ej)]} = eyRsQi^) - hyR&{ey), 
Py = Re{-/[Im(ej.)/ij: -|-Im(/rj)e.v] + [/fvRe(ez) - exR&{hz)]} = hxRe,{e^) - exRe-{h^), 


where Cx, ey, hx, hy are approximately set to be real, and (h^) purely imaginary at Xn = 0 but 
complex at Xn 7 ^ 0. Thus, with the assistance of Eq. (fTTI) . one can easily hnd that Px and Py are 
not 0 for none zero Re(ej.) and Re(/ 2 ^) at Xn 7 ^ 0, which corresponds to the leakage of power 
in the transverse plane. 

It should be noted that such power leakage also leads to difficulties in numerical calculations, 
as seen by the unusual cut-offs in Fig.|2l As Xn increases, the fields of the x-polarized modes 
on the outer-boundary are so large that they cannot be described by the scattering boundary or 
perfectly electric conductor employed in our FEM modal solver. Such a problem can be solved 
by properly including a perfectly matched layer. 


3.4. Perturbation 

It is clear from Eq. (|7]i that the first and second term on the RHS mostly contributes to the 
real and imaginary part variation of j3, respectively. Moreover, the slope of Im(«(.//) against 
X 12 can be calculated by applying perturbation to Eq. (|7]i, assuming that the field prohles at p 
remain the same as those at j3o for tiny change of Axn- This approximation is the essential bit 
in most of the perturbation theory, and is also valid here. Following such logic and substituting 
Hgff = we reach the following result, 

AIm(«,//) jj^^^^eoxhlydxdy 

AZi2 If {eoxhoy - eoyhlJdxdy' 

It is clear from Eq. (fT2l i that the slopes of y-polarization modes are zero because eox ~ 0 and 
liQy Ri 0 SO that fJ^g^^exhQ^,dxdy = 0. Equation (fT2l i shows that the slope of lm{neff) is related 
to the magnitude of how well the in-plane components are confined in the core region. Figure|4] 
shows the calculated slopes by extracting values from Fig. |2lb) (symbols without lines) and 
by our perturbation theory (symbols with lines). It indicates that Eq. (fT^ could follow the 
right trend of simulation results in all cases. The general slope increases with the radius of 
waveguides, since the larger the radius the better the in-plane components are conhned. In 
addition, the slopes obtained from fullwave simulation are slightly larger than our theoretical 
values. The small discrepancy between simulations and theory could be attributed to the fact 
that the contribution of the first term of Eq. Q is neglected in our study. 


3.5. Dispersion of high order modes in bianisotropic waveguides 

Furthermore, to show the impact of bianisotropy on high order modes, the radius of the cylin¬ 
drical waveguide is enlarged to form a multimode waveguide. To keep the discussion simple. 








Fig. 4. Slope of Im(ngyy) of x-polarized modes calculated according to Eq. (112b (symbols 
with lines) as well as from Fig.[2{b) (symbols without lines) versus Xi 2 . Radius of the three 
waveguides are 0 . 12 Ao (black circles), O.I 6 A 0 (blue triangles) and 0.22^ (red diamonds), 
respectively. 


the first high order vector mode group, i.e. TEqi, odd HE 21 , even HE 21 and TMqi, is analyzed. 

Eigure |5la)-13b) show the evolution of Re{neff) and Im{ngff) of the first high order mode 
group as a function of Xn (black dot lines). Eigure |5lc)-13d) shows enlarged mode coupling 
regions in Eig.|3a)-13b). The fundamental modes are shown for reference (blue dashed lines). 
The information in this figure is rich. Eirst, the fundamental modes in this case behave similar to 
that of the single mode case. The mode profiles of the four vector modes at )j(i 2 = 0 as marked by 
el-e4 in Eig.|5lc) are shown in Eig.|5le), corresponding to three different effective mode indices 
because odd HE 21 and even HE 21 modes are degenerate due to circular symmetry. As Xn 
increases, mode coupling occurs between TEqi and odd HE 21 mode as well as even HE 21 mode 
and TMqi mode. This is because linear polarization (in particular y-polarization in the paper) 
is preferred by the bianisotropic medium as Xn increases, the four vector modes evolve into 
linearly polarized modes and the degeneracy between odd HE 21 and even HE 21 mode is broken 
due to required mode coupling. Typical mode profiles of x- and y- polarized odd and even modes 
as marked by fl-f4 in Eig.|5la) are shown in Eig.|5lf), corresponding to X 12 = 0-6 in this case. 
When linearly polarized modes are formed, x-polarized odd and even modes behave similar 
to the fundamental 2 c-polarized modes with components and \\m.{ngff)\ increase, showing 
the modes turn into leaky modes. It is worth noting that the slopes of \m{neff) calculated by 
Eq. (fT^ in this case are close to 1, matching the values extracted from Eig.|5lb). This indicates 
that the in-plane components, i.e., and Hy, are well confined in the core region. As shown 
in Eig.|3c)-l5jd), mode coupling behavior is complicated. However, detailed coupling behavior 
between high order modes is beyond the scope of this paper and are left to future work. 

4. Conclusions 

In closing, we have investigated the mode properties of circularly symmetric waveguides which 
contain bianisotropic media in core layers. Essentially, we find that bianisotropy breaks the po¬ 
larization degeneracy between two linear polarizations, i.e., x- and y-polarizations. The bian¬ 
isotropy turns a mode of undesired polarization (x-polarization in this case) from the degener¬ 
ate pair into leaky mode and leaves the other one unchanged, in contrast to previous findings 
in chirowaveguides or chiral fibers that x- and y-polarization are entangled together. Moreover, 
modes with undesired polarization can be treated as effectively eliminated due to the exponen¬ 
tial attenuation as they propagate. A simple model is built from coupled mode theory to study 
the impact from the bianisotropy, on the mode properties. A perturbation is further applied to 
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Fig. 5. Real part (a) and imaginary part (b) of n^ff versus Xl2- (c) (d) shows enlarged 

region of (a) and (b), respectively, (e) Mode profiles of the first high order mode group 
when X\2 = 0 (f) X\2 = 0.6. Radius of the waveguide is 0.416Ao. 


account for the variation of imaginary part of the effective mode with respect to the chirality 
parameter and shows good agreement with numerical results. It is worthwhile to point out that 
the effective elimination of undesired polarization is general, which also occurs in waveguides 
with non-circular cross-sections, and operates in multimode regime. 

As an outlook, we envisage that the bianisotropic waveguides might be useful in terahertz 
or middle infrared regime. Compact and integrated polarizers is challenging but essential for 
devices and systems operating with a single polarization ll36l438ll . The bianisotropic waveg¬ 
uides proposed here may provide potentials for integrated linear polarizer with high polarization 
selectivity over compact interaction length. Regarding the realizability of bianisotropic waveg¬ 
uide, there are three relevant dimensions: (1) the feature size compatible with current fabrication 
technology, (2) the size of the meta-atoms, (3) the cross-section of waveguide. To respect the 
effective constitutive parameter used here, the size of meta-atoms should be less than I/IOAq. 
Therefore, fabrication of optical bianisotropic waveguides is challenging. However, it may be 
possible to explore direct laser writing (DLW) combined with STED technique ll^ for opera¬ 
tion at middle infrared or longer wavelengths. As the size of SRRs or helix arrays approaches 
half wavelength along the propagation direction, the waveguides would behave similar to ID 
photonic crystals, in which the cross-section of the waveguide shall be large enough to contain 
at least one meta-atom. In this case, we believe it is possible to observe similar effect in the 
pass band. In general, engineering bianisotropic waveguides, or bianisotropic metamaterials, 
may lead to the new possibilities of manipulating light propagation, emission and absorption, 
and interesting phenomenon such as one way transportation of photons. 
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